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Abstract. - The evolution to the steady state of a granular gas subject to simple shear flow is 
analyzed by means of computer simulations. It is found that, regardless of its initial preparation, 
the system reaches (after a transient period lasting a few collisions per particle) a non-Newtonian 
(unsteady) hydrodynamic regime, even at strong dissipation and for states where the time scale 
associated with inelastic cooling is shorter than the one associated with the irreversible fluxes. 
Comparison with a simplified rheological model shows a good agreement. 



O \ Introduction. — A granular gas is a large collection 
i ^o f (mesoscopic or macroscopic) particles which collide in- 
elastically and are usually kept in a state of continuous 
agitation. Apart from their interest in industrial and tech- 
nological applications, granular gases are important at a 
fundamental level as physical systems intrinsically out of 
equilibrium and thus exhibiting a wide spectrum of corn- 
eal plex behavior [1-6] . Although the number of grains in a 
fiuidized granular system is of course much smaller than 



the number of atoms or molecules in a conventional gas, 



■» it is large enough as to make nonequilibrium statistical- 

mechanical concepts and tools applicable. In particular, 
C a kinetic theory approach (based on the Boltzmann and 
I Enskog equations suitably modified to account for inclas- 
^Otic collisions) has proven to be very useful [6]. However, 
^because of the energy dissipation upon collisions and the 
q associated lack of detailed balance and Gibbs equilibrium 
• *state, one is not allowed to take for granted any phe- 
. ^ nomenology that applies to normal gases, unless much 
^ caution is exercised. Quoting Kadanoff, "one might even 
H say that the study of granular materials gives one a chance 



to reinvent statistical mechanics in a new context" [1]. 

Hydrodynamics is one of the key features of standard 
fluids. As is well known, the hydrodynamic description 
of a conventional fluid consists of closing the exact bal- 
ance equations for the densities of the conserved quan- 
tities (mass, momentum and energy) with constitutive 
equations relating the momentum and heat fluxes to the 
conserved densities (usually referred to as hydrodynamic 
fields) and their gradients. If the hydrodynamic gradients 
are weak, the fluxes can be assumed to be linear in those 
gradients, what results in the Navier-Stokes (NS) hydro- 
dynamic description. On the other hand, even if the gradi- 
ents are strong, a (non-Newtonian) hydrodynamic regime 
beyond the NS one is still possible. The conventional 



scenario for the "aging to hydrodynamics" in a normal 
gas can be summarized as follows [7]. Given an arbitrary 
initial state, the evolution proceeds along two successive 
stages. First, during the so-called kinetic stage there is a 
fast relaxation (lasting a few collision times) to a "univer- 
sal" (or "normal") velocity distribution /(r, v;t) that is 
a functional of the hydrodynamic fields (number density, 
flow velocity and temperature). Subsequently, the hydro- 
dynamic stage is described through a slower evolution of 
the hydrodynamic fields as they approach equilibrium or 
an externally imposed nonequilibrium steady state. The 
first stage is sensitive to the initial preparation of the sys- 
tem, while in the hydrodynamic regime the system has 
practically "forgotten" the details of its initial state (ex- 
cept for an implicit dependence on the initial conditions 
through the hydrodynamic fields). If the hydrodynamic 
gradients are small enough when the hydrodynamic state 
is reached, the latter can be described by the NS terms 
in the Chapman-Enskog expansion [8]. However, a nor- 
mal (or hydrodynamic) velocity distribution function is 
not restricted to the NS domain but can apply to the non- 
Newtonian regime as well [9] . 

The applicability of a hydrodynamic description to 
granular fluids is not self-evident at all [1]. In partic- 
ular, the absence of energy conservation gives rise to a 
sink term in the energy balance equation which might 
preclude, except perhaps in quasi-elastic situations, the 
role of the (granular) temperature as a hydrodynamic vari- 
able. At least in the NS regime, however, there exists com- 
pelling evidence from theory [10], simulations [11] and ex- 
periments [12] supporting the validity of a hydrodynamic 
treatment of granular fluids. On the other hand, the in- 
herent lack of scale separation in sheared granular gases 
invalidates a NS description, so that the applicability of 
hydrodynamics in those cases has been controversial [13]. 
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The scenario of aging to hydrodynamics described above 
is known to hold in two limiting cases: finite dissipation 
but no gradients (i.e., in the approach to the homoge- 
neous cooling state) [14] and zero dissipation but large 
gradients [9]. Having that in mind, the question we want 
to address in this paper is, does that scenario still ap- 
ply to granular gases with strong dissipation in the non- 
Newtonian regime? 

Uniform shear flow. — In order to contribute to an 
understanding of the previous question and its possible 
answer, we will focus on inelastic hard spheres in the so- 
called simple or uniform shear flow (USF), which is char- 
acterized by a constant density n, a time-dependent (but 
uniform) granular temperature T(t) and a flow velocity 
with a linear profile u(r) = ayx, a being the constant 
shear rate. There are three basic reasons to choose this 
state to isolate the problem at hand. First, it is macro- 
scopically simple since only a hydrodynamic gradient ex- 
ists (a = du x /dy) and moreover it is a constant. Second, 
this flow possesses a steady state (resulting from the bal- 
ance between inelastic cooling and viscous heating) that is 
inherently non-Newtonian [13,15]. Finally, since both the 
density n and the flow velocity u(r) are independent of 
time, the possible aging to hydrodynamics is enslaved by 
the temporal evolution of the granular temperature, which 
is precisely the quantity casting doubts on the applicabil- 
ity of hydrodynamics. 

In the USF the mass continuity equation is identically 
satisfied, while momentum conservation implies that the 
shear stress P xy is uniform. The energy balance equation 
becomes 



Restricting ourselves to low density granular gases and 
to states that are uniform in the Lagrangian frame of ref- 
erence, so that f(r,v;t) = /(V(r),i) with the peculiar 
velocity V(r) = v u(r), the Boltzmann equation be- 
comes [9, 15] 



(d t - aV y dvJ f(V,t) = J[V\f(t)}, 



(2) 



where J[V|/(t)] is the (inelastic) Boltzmann collision op- 
erator. Equation (2) must be complemented with an 
initial condition /(V, 0) = /°(V), so that the solution 
is actually a functional of the initial distribution, i.e., 
f(V,t) = f(V,t\f°). Analogously, the velocity moments 
of /, such as the pressure tensor Pij(t\f°) are also func- 
tional of f°. Since the only time-dependent hydrody- 
namic variable is the temperature T(t) and the only hy- 
drodynamic gradient is the shear rate a, the existence of a 
hydrodynamic regime implies that, after a certain number 
of collisions per particle, 

/(V, t\f) n [m/2T(t)f 2 /* (C(i); a*(t)), (3) 



8 t T(t) = -(2a/3n)P xy (t) - ((t)T(t) 



(1) 



where ((t) is the cooling rate due to inelasticity. Given 
a shear rate a and a coefficient of normal restitution a, 
a stationary temperature T s = lim^oc T(t) is reached 
when the viscous heating and the inelastic cooling can- 
cel each other. This steady state can be conveniently 
described by introducing dimensionlcss quantities. Let 
u(T) = nT/rj a (T) oc nT 1 / 2 , where r] (T) is the NS shear 
viscosity in the elastic case, denote an effective collision 
frequency. We can then define a reduced shear rate a* (t) = 
a/v(T(t)), a reduced cooling rate £*(*) = ((t) /v(T(t)) and 
a reduced shear stress P* (t) = P xy (t) / nT (t) . Thus, one V*( a *) = - 



where C(t) = V/ y / 2T(t)/m is the (peculiar) velocity in 
units of the (time-dependent) thermal speed, m being the 
mass of a grain. The scaled velocity distribution func- 
tion f*(C;a*) is, for a given value of the coefficient of 
restitution a, a universal function in the sense that it is 
independent of the initial state / and depends on the 
applied shear rate a through the reduced scaled quantity 
a* only. In other words, if a hydrodynamic description 
is possible, the different solutions /(V, t\f°) of the Boltz- 
mann equation (2) would be "attracted" by the universal 
form (3). For asymptotically long times, the steady state 
would eventually be reached, i.e., f*(C;a*) —> /*(C) = 
/*(C; a*). Equation (3) has its counterpart at the level of 
the velocity moments. In particular, the pressure tensor 
Pij(t\f°) would become 



P«(t|/°)^nT(i)i*(a*(t)) 



(4) 



with universal functions P*j(a*). The non-Newtonian 
character of the hydrodynamic regime can be character- 
ized by the nonlinear (reduced) shear viscosity rj*(a*) and 
viscometric function 'J* (a*) defined by 



has §a*|P* y J = (* in the steady state. The stationary 
temperature is approached from below (T(t) < T s ) if the 
shear rate is large enough and/or the inelasticity is small 
enough so that the heating rate (due to viscous effects) 
^a* (t)\P* y (t)\ prevails over the cooling rate (due to inelas- 
ticity) C*(^)i we w iU refer to this situation as the heating 
case. Otherwise, the stationary temperature is reached 
from above and this will be referred to as the cooling case. 
In this context, the question posed above can be rephrased 
as, does an unsteady hydrodynamic regime exist in both 
the heating and cooling cases before the steady state is 
achieved? 



a* 



#*(a*) - 



P* xx {a*) - P; y (a* 

„*2 



(5) 



Rheological model. — Before presenting the simula- 
tion results, it is worth considering a simple kinetic model 
[16] in which the Boltzmann collision operator J[V|/(t)] 
is replaced by a relaxation-time term toward the local 
equilibrium distribution plus a drag-force term mimick- 
ing the cooling effects due to inelastic collisions. Applied 
to the USF, it closes eq. (1) with the evolution equa- 
tions [15,17,18]: 

d t P xy = -aP yy - f3v(T)P xy - C(T)P xy , (6) 

d t P yy = -0v(T) [Py y - nT) - C(T)P yy , (7) 
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where C(T) = CM T ) with C = M 1 ~ " 2 ) and hcre we 
take (3 = |(1 + a). Equations (6) and (7) can also be 
obtained, with (3 = ^(1 + a) (2 + a), from the Boltzmann 
equation in Grad's approximation [15, 19]. 

The set of coupled equations (1), (6) and (7) must be 
solved numerically [15], the solution including both the 
transient kinetic regime and the hydrodynamic evolution 
stage. In order to extract the hydrodynamic solution in 
an analytical way, we will introduce an additional simplifi- 
cation similar to that carried out in the elastic case [9,20]. 
As mentioned before, v(T) oc nT q with q = \. How- 
ever, here we will temporarily view 5 as a (small) free 
parameter, so that the solutions to eqs. (1), (6) and (7) 
depend parametrically on q [21]. If q — 0, then v = const 
and a* — const, so there is no steady state (except if 
a* takes a specific value). However, for sufficiently long 
times the scaled quantities P*j reach well defined (hydro- 
dynamic) values which depend on a* and not on the initial 
state [22]. Next, carrying out a first-order perturbation 
analysis around q = and constructing Pade approxi- 
mants, one finally gets [22] 

rj*(a*) = /3- 1 [l + 2 7 (a*)]- 2 {l + i[C7/5-2 7 (a*)] 
x[l-6 7 (a*)]/[l + 6 7 (a*)] 2 }~\ (8) 

L 2[l + 6 7 (a*)J > 

r C7/3 - 2 7 (a*) 1 - 1 
x i 1 wm n_M g 

I [l + 6 7 (a*)] 2 J 

where 7 (a*) = § sinh 2 [± cosh _1 (l + 9a* 2 //? 2 )] and we 
have already set q — \. Equations (8) and (9) provide 
the hydrodynamic non-Newtonian viscosity and viscomet- 
ric function as explicit functions of the reduced shear rate 
a* (and of the coefficient of restitution a) in our simplified 
rheological model. The steady state corresponds to just 
one point on the curves n*{a*) and 'J* (a*). In the model, 
it is given by a* = ^/iC*/W{P + C) [15,17,18], which ver- 
ifies the condition 7 (a*) = (*/20, and so 77* = i]*(a* s ) = 
P/(P + C) 2 and = **(a*) = 2/3/((3 + C*) 3 . If we 
formally take the limit of vanishing shear rate, we recover 
the NS viscosity of the inelastic gas predicted by the model 
[18], namely lim a .^ V*( a *) = {0 + C*/ 2 )~\ as well as the 
Burnett value lim o .^ = 2/(/3 + C*/2)(/3 + C*). We 

have checked that, for a > 0.5, Eqs. (8) and (9) yield val- 
ues which are hardly distinguishable from those obtained 
by a numerical solution of the set (1), (6) and (7) [22]. 

Monte Carlo simulations. — In order to test 
whether the hydrodynamic regime (3) exists or not, we 
have solved the Boltzmann equation (2) by means of the 
direct simulation Monte Carlo (DSMC) method [23]. As in 
ref. [24], we have taken N = 10 4 simulated particles and a 
varying time step St = 10" 3 r° ^JWJT, where T° is the ini- 
tial temperature and r° = A/^2T°/m ~ QMttv- 1 ^ ) 



is the initial mean free time, A being the mean free path. 
In order to improve the statistics, the results are averaged 
over 100 independent realizations. The values of the co- 
efficient of restitution considered have been a = 0.5, 0.7 
and 0.9. For each value of a, four different shear rates 
have been taken: a = 0.01/r°, a = 0.1/r°, a = 4/r° 
and a = 10/r°. The two first values of the shear rate 
(a = 0.01/r and a = 0.1/r°) are small enough to cor- 
respond to cooling cases, even for the least inelastic sys- 
tem (a = 0.9), while the other two values (a = 4/t° and 
a = 10/r ) are large enough to correspond to heating 
cases, even for the most inelastic system (a = 0.5). For 
each one of the twelve pairs (a, a), five widely different 
initial conditions have been chosen, so that sixty indepen- 
dent states have been simulated. The first initial condition 
(here referred to as A) is just the (local) equilibrium one, 
namely 

/°(V) =n(m/27rT°) 3/2 e- my2 / 2T °. (10) 

The other four initial conditions are of the anisotropic 
form 

/°(V) = (n/2)(m/27rT°) 1/2 e- mV ?/ 2T ° 

x [5 (V x - V° cos <j>) 8 (V y + V° sin 0) 
+5(V x + V°cos(f))S(V y - V°sin0)] (11) 

where S(x) is Dirac's distribution, V° = ^/2T°/m is the 
initial thermal speed and 4> G [0, ir] is an angle charac- 
terizing each specific condition. The pressure tensor cor- 
responding to eq. (11) is given by P xx = 2nT° cos 2 4>, 
P° y = 2nT°sin 2 </), P° 2 = nT° and P° y = -nT°sin20. 
The four values of (f> considered are (f> — kir/A with k = 0, 
1, 2 and 3; we will denote the respective initial conditions 
of type (11) as B0, Bl, B2 and B3. Note that P° y = in 
the initial conditions A, B0 and B2, while P xy = — nT° in 
the case of Bl. On the other hand, P xy takes an "artifi- 
cial" positive value (P xy = nT°) if the system is initially 
prepared with the condition B3. According to eq. (1), this 
implies that during the very first times the fluid undergoes 
a double cooling effect: the inelastic cooling plus an ar- 
tificial shear stress cooling. It takes some time before a 
"natural" negative value of P X y(t) is established and the 
usual viscous heating effect competes against the inelastic 
cooling. 

Results and discussion. — The solution of the 
Boltzmann equation (2) through the DSMC method al- 
lows one to monitor the temporal evolution of velocity 
moments, such as the temperature T(t) and the pressure 
tensor Pij(t), as well as of the velocity distribution func- 
tion /(V, t) itself. In particular, one can follow the re- 
duced shear rate a*(t) = a/u(T{t)) (which decreases in 
the heating states and increases in the cooling states) and 
the reduced shear viscosity n*(t) = —P xy (t)/nT(t)a*(t). 
A parametric plot of 77* (t) vs a* (t) is then useful to test 
the establishment of a hydrodynamic regime in which 
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Fig. 1: (Color online) Reduced shear viscosity rf{t) (top panel) 
and number of collisions per particle (bottom panel) versus 
the reduced shear rate a*(t) for a = 0.5 in the two heating 
cases a = 4/r° (red lines) and a = 10/r° (blue lines). In the 
top panel, the circle represents the steady-state point (oj,^), 
while the thin solid line corresponds to the hydrodynamic func- 
tion, eq. (8), obtained from our simplified rheological model. 
Note the logarithmic scale in the vertical axis of the bottom 
panel. 



77* (t) — > r]*(a*(t)), where the function rj*(a*) must be in- 
dependent of the initial conditions. Such a parametric plot 
is shown in the top panel of fig. 1 for the most inelastic 
gas (a = 0.5) and for the two heating cases (a = 4/r° 
and a — 10/r°). The curve representing the rheologi- 
cal model (8) is also included. It is quite apparent that 
the ten curves are attracted to a common universal curve 
which, in addition, turns out to be excellently described 
by our simplified model. The bottom panel of fig. 1 dis- 
plays the temporal evolution of a*(t), as measured by the 
accumulated number of collisions per particle (i.e., the to- 
tal number of collisions in the simulations divided by the 
total number of particles) . We can observe that the states 
with the highest heating effect (a = 10/t°) reach the hy- 
drodynamic stage after about 1 collision per particle only, 
while this aging period takes a little longer (about 2 col- 
lisions per particle) in the states with a = 4/r°. Once 
the state joins the hydrodynamic curve, it moves along it 
until the steady state is reached, what occurs after typi- 
cally 10 collisions per particle since the initial state. This 
means that about 80-90% of the duration (measured by 
the number of collisions) of the evolution to the steady 
state is occupied by the hydrodynamic stage. It is also in- 
teresting to remark that, for a given value of a, the slowest 
heating takes place for the initial condition B0, followed 
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Fig. 2: (Color online) Same as in fig. 1, but for the two cooling 
cases a = 0.01/r° (blue lines) and a = 0.1/r° (red lines). Note 
the logarithmic scale in the horizontal axes and the linear scale 
in the vertical axis of the bottom panel. 



by B3, A and Bl, while the fastest heating corresponds 
to B2. This implies that the state starting from B0 is the 
one joining the hydrodynamic curve at a larger value of a* 
(a* ~ 7 and a* ~ 3 for a — 10/r° and a — 4/r°, respec- 
tively), while the one starting from B2 does it at a smaller 
value (a* ~ 3 and a* ~ 1.5 for a = 10/r° and a = 4/r°, 
respectively). 

Let us consider now the cases a = 0.01/t° and a = 
0.1/r°. The temporal evolution is dominated now by the 
inelastic cooling (i.e., the time scale C _1 associated with 
the cooling is shorter than the time scale \nT j a\P xy \ asso- 
ciated with the viscous heating) and so the existence of a 
hydrodynamic regime might appear as more doubtful than 
in the heating cases. That this is not actually the case can 
be concluded from fig. 2. Again, the ten curves tend to 
collapse to a common curve and again the latter practi- 
cally coincides with the theoretical prediction (8), except 
perhaps near the NS region of small a*. In contrast to 
the heating cases depicted in fig. 1, however, the tempo- 
ral evolution of T(t), and hence of a*(t) oc T~ 1//2 (t), is 
practically independent of the initial condition, especially 
in the case a = 0.01/r°. This is a consequence of the 
fact that for these low values of ar° the viscous term in 
eq. (1) can be neglected for short times and so the tem- 
perature initially evolves as in the homogeneous cooling 
state (decaying practically exponentially with the number 
of collisions), hardly affected by the details of the initial 
state. By the time the viscous heating effect becomes com- 
parable to the inelastic cooling, the hydrodynamic regime 
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Fig. 3: (Color online) Reduced shear viscosity 77* (a*) (top 
panel) and viscometric function <3/* (a*) (bottom panel) for the 
hydrodynamic part of the evolution to the steady state (repre- 
sented by a circle) for, from top to bottom, a = 0.5, 0.7 and 
0.9. The color code is the same as in figs. 1 and 2. The thin 
solid lines correspond to the hydrodynamic functions, eqs. (8) 
and (9), obtained from our simplified rheological model. 

has already been attained. It is important to note that 
the number of collisions needed to reach the hydrody- 
namic behavior (about 5 and 10 collisions per particle 
for a = 0.01/t and a = 0.1/r°, respectively) is longer 
in the cooling cases than in the heating cases. On the 
other hand, the total duration of the evolution period is 
also longer than in the heating cases and so the granular 
gas still lies on the hydrodynamic curve for most of the 
evolution time (as measured by the number of collisions 
per particle). The duration of the transient kinetic stage 
in the cooling cases (about 5-10 collisions per particle) is 
consistent with the one observed in molecular dynamics 
simulations of the homogeneous cooling state [14]. 

We have observed behaviors similar to those of figs. 1 
and 2 for the viscometric function and for the other in- 
elasticities (a — 0.7 and a — 0.9). While the number of 
collisions the system needs to loss memory of its initial 
state is practically independent of a, the total duration of 
the evolution period tends to increase with a, especially 
in the cooling cases [22] . This means that the more elastic 
the system, the smaller the fraction of time spent in the 
kinetic regime prior to the hydrodynamic stage. 

The common hydrodynamic parts of the functions 
r]*(a*) and $*(a*) for the ten heating cases and the ten 
cooling cases are displayed in fig. 3 for each one of the 
three inelasticities considered. This corresponds to the 




Fig. 4: (Color online) Marginal distribution functions g x (C x ) 
and g y (C y ) for a = 0.5 and a* = 0.40, 0.66 and 1.26. The color 
code is the same as in figs. 1 and 2. 

shear rate windows 0.4 < a* < 1.3, 0.3 < a* < 1.3 and 
0.2 < a* < 1.2 for a = 0.5, 0.7 and 0.9, respectively. The 
degree of overlapping of the curves is excellent, although 
the statistical fluctuations inherent to the DSMC method 
are more important for ^*(a*) than for 77* (a*), especially 
in the cooling branch of the curves. It is also apparent that 
our simplified model, eqs. (8) and (9), describes reasonably 
well the rheological properties of the sheared granular gas 
Figures 1-3 confirm eq. (4), i.e., the existence of well 
defined hydrodynamic rheological functions P*j(a*) act- 
ing as attractors in the evolution of the pressure tensor 
Pij(t\f°), regardless of the initial preparation /°. This in 
turn provides indirect support to the stronger statement 
(3). To test it in a more direct way, we have considered 
the marginal distributions [24] g x (C x ;a*) and g y (C v ] a*), 
defined by 

/oo />oo 
dC z / dC VtX f*(C;a*), (12) 
-00 Jo 

for a — 0.5 and three values of a* that, according to figs. 1- 
3, are within the hydrodynamic range, namely a* = 0.40, 
0.66 (cooling branch) and 1.26 (heating branch). The re- 
sults are plotted in fig. 4, which shows an excellent over- 
lapping of the ten curves for each value of a*, although 
obviously the tails of the distribution exhibit larger statis- 
tical fluctuations than the thermal region. As expected, 
the anisotropic features of the velocity distribution in- 
crease with the shear rate. The steady state corresponds 
to a* — 0.92 [24] and hence the shapes of the associated 
steady distributions (not shown) lie in between those for 
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a* = 0.66 and a* = 1.26. 

In conclusion, the results reported in this paper, along 
with a more extensive study that will be published else- 
where [22], strongly support that the conventional sce- 
nario of aging to hydrodynamics remains essentially valid 
for granular gases, even at high dissipation, for non- 
Newtonian states and for situations where the time scale 
associated with inelastic cooling is shorter than the one 
associated with the irreversible fluxes. 
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